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^NJ We present a stochastic method for the calculation of baryon three-point 

O functions that is more versatile compared to the typically used sequential 

method. We analyze the scaling of the error of the stochastically evaluated 
three-point function with the lattice volume and find a favorable signal-to- 

• ^ noise ratio suggesting that our stochastic method can be used efficiently at 

^ large volumes to compute hadronic matrix elements. 
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1 Introduction 



Hadron structure calculations in lattice QCD have emerged as a powerful tool for pro- 
viding comparison to and guidance for experiments, see e.g. refs. [H El [3l H] . Examples 
include moments of generalized parton distribution functions, as well as form factors. 
Lattice computations of such quantities have been carried out at several values of the 
lattice spacing allowing the continuum limit to be taken. In addition, small, sometimes 
even physical, values of the pion mass have been employed in the calculation of these 
quantities leading to an improved understanding of their quark mass dependence and 
how they approach the physical point. Unfortunately, studies of excited state contribu- 
tions O El El El [9] suggest that for some quantities these effects can play an important 
role as a systematic uncertainty that affects hadronic three-point function computations. 
Safely accounting for these effects requires large statistics, hence methods to speed-up 
these calculations are highly desirable. 

The progress in nucleon matrix element calculations on the lattice has prompted an 
effort to go beyond the simplest observables and to pursue a larger variety of interesting 
hadronic quantities. The evaluation of more observables will deepen our knowledge of 
hadron structure and provide a more comprehensive test of QCD. However, using the 
conventional sequential method [10] to calculate these matrix elements, it is necessary 
to perform a new computation of the needed quark propagators for each observable of 
intereslj^ This then leads to a high computational demand if many physical quantities 
are being sought. 

In this paper, we describe an alternative approach, based on a stochastic method, 
that allows us to obtain a large class of observables with only a single computation of 
the propagators. To this end, we employ stochastically computed all-to-all propagators. 
Since a calculation based on a stochastic evaluation of propagators may lead to very 
noisy results, we perform a detailed study to determine whether the stochastic noise can 
be controlled with a moderate number of stochastic sources. We determine the signal- 
to-noise ratio as a function of the lattice size to test whether our stochastic method 
can be used in large volumes, such as 48'^ x 96, that are used in contemporary lattice 
computations. These questions are addressed specifically for the example of the axial 
charge of the nucleon. We would like to emphasize that we do not perform an analysis 
of systematic effects, since our goal is solely to test the stochastic method. 

During the course of our work, a similar stochastic approach was employed in the 
calculation of meson three-point functions [11], where a (heavy) all-to-all propagator is 
estimated stochastically. There it was found that the stochastic method is competitive 
and in some cases even superior to the sequential one. Of course, it is not guaranteed 
that the same conclusions hold for baryon matrix elements, since those are subject to a 
stronger exponentially decreasing signal-to-noise ratio [12j. 

This paper is organized as follows: in Sec. [2] we outline our stochastic method, in Sec. [3] 

^The sequential method allows to compute either the matrix element of a single hadron or the matrix 
elements of a single current with one sequential inversion. In particular, a new computation of 
propagators has to be performed, when computing matrix elements of another hadron (fixed sink 
method) or another current (fixed current method). 
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we present the results of this method and compare to those obtained by the fixed sink 
method, and we summarize our findings in Sec. |4j 



2 Stochastic method for baryon three-point functions 

Quantities that are needed for investigating hadron structure can be extracted by com- 
puting matrix elements of baryons with local operators. In lattice QCD these baryon 
matrix elements are obtained from baryon three-point functions in Euclidean space-time 
that are of the form 



x,y 



e-^P-y{B{x,t)\0{y,T)\m) , (1) 
0{y,T)=q{y,T)Tq{y,T) . (2) 



r represents a combination of 7-matrices and covariant derivatives and we have used 
translational invariance to set the source point to zero. Naively one would need an all- 
to-all propagator from every lattice point (y, r) to all points (x, t) for the evaluation of the 
above three-point function, which is, of course, prohibitively expensive to calculate. Such 
a demanding computation can be circumvented by applying the sequential method to 
perform the summation over the spatial coordinates of either the sink or the current [lOj . 
For the example of the fixed sink method, the momentum as well as the time slice of 
the sink are fixed and an additional inversion for each flavour is needed. An alternative 
approach is to estimate the all-to-all propagator stochastically, which is the method that 
we explore in this paper. 

A generic three-point function of a baryon B is defined as 

Ib is the baryon interpolating field and A and A' summarize the indices depending on 
the quantum numbers of the baryon B, which are appropriately contracted with the 

( B) 

function Cyiy/'- '^^^ insertion time of the operator is denoted by r. For illustration let 
us now consider the three-point function of a proton and the operator dTd. We use the 
interpolating field Ip possessing the quantum numbers of the proton, namely 

where C = ^7072 is the charge conjugation operator. In terms of quark propagators, the 
connected piece of this three-point function reads 



Cj^) (t,r;p,p-';P) =P""'^ e-^P-ye-^^'^ e,ta Sa'b'c' {Cj,)^,y {C-f,);^, T^'^d^'d 

M „.\ aid) 
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where the up (down) quark propagator is denoted by S*^"^ {S^'^^). P is an appropriate 
spin projector, which we will specify later. 

The sequential method with fixed sink makes use of the fact that we can perform the 
sum over x in Eq. ([T]) by an additional inversion. Then a generalized propagator for fixed 
time slice, projector P and sink momentum f/ is obtained, as indicated by the shaded 
area in the left diagram of Fig. [TJ This renders the explicit calculation of an all-to-all 
propagator unnecessary. 

Our alternative method uses a stochastic estimate of the all-to-all propagator appear- 
ing in three-point functions like Eq. ([s]). Such an estimate is obtained via 



,x) rir{x) 



where Ai is the Dirac matrix. In the above equations we have suppressed Dirac and 
color indices, rjr is a random source obeying 

Ns 



Ns 



r=l 



The stochastic method is diagrammatically illustrated in the right diagram of Fig. [T] 
Using the stochastic estimate, we can decompose the double sum in Eq. ^ into the 
product of two single sums, which is significantly computationally cheaper than the 
naive double sum and reads 

(t,r;p,p';P) =P""X ^"'''''^cbaea'b'c' (Cts)^^ {Cj,\fs r^'^'^d'd 

X 

y 
y 

We have suppressed the average over the number of stochastic samples, cf. Eq. Q, and 
used the 75 Hermiticity property of the Dirac matrix to obtain the above expression. As 
before, we use superscripts to denote the quark flavor. 

The drawback of the stochastic method is that we have to average over a sample 
of Ns stochastic sources. This requires Ns inversions compared to just twelve (one 
inversion per Dirac and color index) in the sequential method. However, a major benefit 
of this method is its flexibility, since we do not need to fix the spin projector or the sink 
momentum, nor even the sink time slice, in principle. 
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Figure 1: Diagrammatic illustration of the sequential method through the sink (left) 
and the stochastic method (right). 

3 Assessment of the stochastic method 



To test the applicability of the stochastic method, we need to determine how large Ns 
must be in order to keep the stochastic noise under control. This depends on the observ- 
able of interest. To be concrete, we compute a relatively simple benchmark observable 
of nucleon structure, namely the nucleon axial charge qa, using our stochastic method. 
This quantity can be obtained from a nucleon matrix element of the isovector axial cur- 
rent. For the evaluation of matrix elements of the nucleon, we need to introduce the 
zero momentum nucleon two-point function, 



cf)(t) = 5:Tr 



^(l+7o)(xW(^',0 ^^^^(0) 



We then examine the ratio Rg^ of the nucleon three-point and two-point correlation 
functions. 



cr\t) 



(5) 



c'i!Xii,r;p = i^ = Q) = Y.Tr rfe(x(^)(x,t) Ofc(y,r)xW(0) 



Ok (y, t) = u{y, T)j5-fku{y, t) - d{y, r)757fc(i(y, r), 
i 

Tfc = ;^ (1 + 7o)757A:, A; =1,2,3 



In the limit of large Euclidean time separations, Rg^ converges to qa up to the renor- 
malization factor Za, 



ZARgA{'t,T 



gA for t — )• oo, T — )• oo and (t — r) 



oo. 



In this paper we use the value of the renormalization constant Za = 0.757(3) determined 
non-perturbatively jlSl E] . 

In order to demonstrate that the stochastic method can indeed produce results with 
a reasonable computational effort and is potentially competitive with the sequential 
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method, we performed a benchmark calculation with Nf = 2 + 1 + 1 flavors of quarks. We 
employed twisted mass fermions at maximal twist with a lattice spacing of a ~ 0.086 fm 
determined from the nucleon mass, a pion mass of 111^^ ~ 380 MeV and a volume of 
~ (2.7 fm)^. In Fig. [2j we show Rg^{t,T) obtained using the stochastic method as a 
function of the insertion time r for a fixed source-sink separation t = 12a. We compare 
to the value Rgj^{t = 12a, t = 6a) obtained using the sequential method. For different 
values of r close to the middle of the plateau the picture is similar. We use spin-color 
diluted random 2'(4) vectors as stochastic sources, 

<(f , t) = 5a,aoSa,aoSt,to^ix), Oq G {0, 1, 2}, Oq G {0, 1, 2, 3}, fj{x) £ Z(4) . 



We have used a fixed number of gauge field configurations N, 
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Figure 2: Plateau region of the ratio Rgj^{t = 12a, r) obtained from the stochastic 
method on one Nf = 2 + 1 + 1 ensemble with a pion mass of m-Tr ~ 380 MeV 
and a lattice spacing a « 0.086 fm. We use A''^ = 2, 4 and 6 spin-color diluted 
stochastic time-slice sources, with the source located at the sink. The source- 
sink separation is 12a and we compare to the standard sequential method, of 
which we show the ratio Rgj^{t = 12a, t = 6a) with the light gray band. The 
dark gray band indicates the PDG value 1 15]. 

stochastic and sequential approaches. Our observation is that using at most Ns = 4 
spin-color diluted stochastic noise vectors per configuration for the estimate of the all- 
to-all propagator is sufficient to reach the same statistical accuracy as with the sequential 
method. 

In terms of inversions, Ns = 1 corresponds to using the same number of inversions as 
in the sequential method i.e. (4 x 12) per gauge field configuration, for every additional 
set of stochastic sources we need (2 x 12) additional inversions to obtain the forward 
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and back propagators, respectively. Thus this method would require four times more 
inversions. We would like to remark, however, that in the stochastic approach we can 
compute the correlation functions of proton and neutron without additional inversions, 
in contrast to the sequential method. In addition we used 3 operators k = 1,2,3 in 
Eq. ([5]), and correspondingly 3 spin projectors for the stochastic method, again without 
the need of additional inversions when using the stochastic method. Our observation is 
that this procedure reduces the statistical error by about a factor \/6, corresponding to 
a factor of about 6 in statistics. We would like to remark however that this effective 
gain in statistics is rather specific for gA and will change for other observables. 

Having demonstrated that it is in fact possible to compute gA using the stochastic 
method with a reasonable computational effort compared to the sequential method, we 
would like to know how the situation changes when the volume is varied. A potential 
danger of our stochastic method is that the number of stochastic sources required to 
reach the same precision as the sequential method may increase rapidly as the volume 
increases. 

To study the volume effects, we use Nf = 2 flavors of maximally twisted mass fermions, 
instead of the A'^^ = 2 + 1 + 1. We expect that the stochastic noise should not noticeably 
depend on the number of dynamical flavors, and for the Nf = 2 case, there exists a 
series of four different volumes at the same value of the lattice spacing, a ~ 0.089 fm |16] 
and a pion mass, 171-,^ ~ 300 MeV |l7j. These volumes are V = x T, where L/a = 
16, 20, 24, 32 and the temporal extent of the lattice is T = 2L in all cases. This enables 
us to thoroughly study the volume dependence of the stochastic method over a relatively 
large range of spatial volumes, from about (1.4 fm)^ to (2.8 fm)^. This corresponds to 
1.95 < m^L < 3.9. 

We performed an analysis using the stochastic method on a flxed number of gauge 
fleld conflgurations A^gauge = 200 at each of the four volumes. The source-sink separation 
is fixed to 12a in all cases, which corresponds to about 1.067 fm. 

In Fig. [3] we show the relative error of gA as a function of the number of stochastic 
sources for two of the four volumes L/a = 24 and L/a = 32, where also the sequential 
method has been applied. For both volumes a convergence towards the error of the 
sequential method can be seen, which looks better for the larger volume, where the error 
is close to the error of the sequential method for Ns = 4, to be conservative. This 
confirms the observation in the Nf = 2 + 1 + 1 calculation mentioned above, which was 
done at about the same physical volume. The error of the error is not shown, but it is 
roughly of order 10% of the error. For smaller volumes of the lattice, the conversion is 
worse, but note that statistical fiuctuations are larger in this case. 

We show the combined statistical and stochastic error oi gA obtained using the stochas- 
tic method with A^5 = 8 for four different volumes in Fig. [4j In order not to be subject to 
a systematic error we do not fit the ratio given in Eq. ((5j) using an estimated plateau range 
r between source and sink. Instead, we take the renormalized ratio Rgjy{t = 12a, r = 6a) 
and its error in the middle between source and sink as our estimate for gA, where con- 
tributions from excited states are expected to be the smallest. For the larger volumes 
the error scales like the one of the sequential method, V~^'^ [18]. Therefore the plot is 
consistent with the error being dominated by the gauge noise for larger volumes, which 
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Figure 3: The filled circles show the relative error on the value obtained for as a 
function of the number of spin-color diluted stochastic sources Ns per gauge 
field configuration. In the left panel we show results for a lattice extent of 
L/a = 24 and in the right panel we show results for L/a = 32. The dashed 
lines represent the error when using the sequential method. A fixed number 
of gauge field configurations A'gauge = 200 was used, with a pion mass of 
m,r ~ 300 MeV and a lattice spacing of a « 0.089 fm. 



we have demonstrated above, see Fig. [3| Thus this method appears to work as well or 
even better at the larger volumes typically employed in current calculations, an obser- 
vation which, of course, needs to be verified for other quantities and different physical 
situations. 



4 Conclusion 

We have applied a stochastic method for the calculation of nucleon matrix elements using 
spin-color diluted time slice sources. We have taken the case of the nucleon axial charge 
as a typical example of a three-point function to explore the method. Our conclusions for 
our test case is that the error is comparable to the error of the sequential method already 
at a moderate number of stochastic sources, namely four spin and color diluted timeslice 
stochastic sources, when using the same number of gauge field configurations. In this 
particular case we have effectively increased the statistics used in the stochastic method 
by a factor 6 through averaging over neutron and proton correlators and using 3 different 
currents, which does not require the computation of new propagators. Moreover, our 
results indicate that the convergence behaviour in the number of stochastic sources, Ns, 
appears to improve when the volume is increased. 

Since the stochastic method needs of 0(10) more inversions (we needed Ns = 4 for 
but this maybe different for other matrix elements) it is competitive with the sequential 
method when one computes 0(10) more types of matrix elements. In the case of g^ 
the increase in computational effort is easily compensated by computing 6 different 
matrix elements with 3 different spin projections for the proton and the neutron. Thus, 
even with the additional computational overhead, the great versatility of the stochastic 
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Figure 4: The filled circles show the error of qa using the stochastic method with a 
fixed number of stochastic sources Ns = 8 on a fixed number of gauge fields 
-/Vgaugc = 200 for four different lattice sizes. In all calculations, the pion mass 
is TOtt ~ 300 MeV and the lattice spacing is a w 0.089 fm. The dashed lines 
indicate a scaling proportional to and and are solely meant to guide 

the eye. The error of the error is not much bigger than the symbol size, hence 
we do not show it. Please note the double logarithmic scale. 

method explained in this paper outweighs the sequential method when many baryon 
matrix elements or form factors are computed. 
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